Loading required package: Rcpp
Loading 'brms' package (version 2.20.4). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: 'brms'
The following object is masked from 'package:stats':
ar
library(modelr)
Attaching package: 'modelr'
The following object is masked from 'package:broom':
bootstrap
library(tidybayes)
Attaching package: 'tidybayes'
The following objects are masked from 'package:brms':
dstudent_t, pstudent_t, qstudent_t, rstudent_t
library(ggridges)
Attaching package: 'ggridges'
The following objects are masked from 'package:tidybayes':
scale_point_color_continuous, scale_point_color_discrete,
scale_point_colour_continuous, scale_point_colour_discrete,
scale_point_fill_continuous, scale_point_fill_discrete,
scale_point_size_continuous
# devtools::install_github("seananderson/ggsidekick") ## not on CRAN library(ggsidekick); theme_set(theme_sleek())# Load functionshome <- here::here()fxn <-list.files(paste0(home, "/R/functions"))invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))
d <-#read_csv(paste0(home, "/data/for-analysis/dat.csv")) %>% read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv") %>%filter(age_ring =="Y", # use only length-at-age by filtering on age_ring!area %in%c("VN", "TH"))
Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 43,886 rows (13%), 294,574 rows remaining
# Minimum number of observations per area and cohortd$area_cohort <-as.factor(paste(d$area, d$cohort))d <- d %>%group_by(area_cohort) %>%filter(n() >100)
# Minimum number of observations per area, cohort, and aged$area_cohort_age <-as.factor(paste(d$area, d$cohort, d$age_bc))d <- d %>%group_by(area_cohort_age) %>%filter(n() >10)
# We can also consider removing individuals where the SE of k is larger than the fitIVBG %>%#mutate(keep = ifelse(k > quantile(k_se, probs = 0.95), "N", "Y")) %>%mutate(keep =ifelse(k < k_se, "N", "Y")) %>%#filter(row_id < 10000) %>%ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +geom_point(alpha =0.5) +facet_wrap(~keep, ncol =1) +geom_errorbar(alpha =0.5) +NULL
mutate: new variable 'keep' (character) with 2 unique values and 0% NA
# Add back cohort and area variablesIVBG <- IVBG %>%left_join(d %>%select(ID, area_cohort)) %>%separate(area_cohort, into =c("area", "cohort"), remove =TRUE, sep =" ") %>%mutate(cohort =as.numeric(cohort))
Joining with `by = join_by(ID)`
left_join: added 2 columns (area_cohort_age, area_cohort)
> rows only in x 0
> rows only in y (146,393)
> matched rows 142,023 (includes duplicates)
> =========
> rows total 142,023
mutate: converted 'cohort' from character to double (0 new NA)
# Compare how the means and quantiles differ depending on this filteringIVBG_filter <- IVBG %>%drop_na(k_se) %>%#filter(k_se < quantile(k_se, probs = 0.95)) %>% filter(k_se < k)
group_by: 2 grouping variables (cohort, area)
summarize: now 306 rows and 6 columns, one group variable remaining (cohort)
ungroup: no grouping variables
# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.
Rows: 380 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): area, model
dbl (2): year, temp
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (cohort)
VBG <- VBG %>%left_join(pred_temp, by =c("area", "cohort"))
left_join: added 2 columns (temp, model)
> rows only in x 0
> rows only in y ( 74)
> matched rows 306
> =====
> rows total 306
# Save data for map-plotcohort_sample_size <- IVBG %>%group_by(area, cohort) %>%summarise(n =n()) # individuals, not samples!
group_by: 2 grouping variables (area, cohort)
summarise: now 306 rows and 3 columns, one group variable remaining (area)
VBG <-left_join(VBG, cohort_sample_size, by =c("cohort", "area"))
left_join: added one column (n)
> rows only in x 0
> rows only in y ( 0)
> matched rows 306
> =====
> rows total 306
write_csv(VBG, paste0(home, "/output/vbg.csv"))# Calculate the plotting order (also used for map plot)order <- VBG %>%ungroup() %>%group_by(area) %>%summarise(min_temp =min(temp)) %>%arrange(desc(min_temp))
ungroup: no grouping variables
group_by: one grouping variable (area)
summarise: now 10 rows and 2 columns, ungrouped
order
# A tibble: 10 × 2
area min_temp
<chr> <dbl>
1 SI_HA 7.83
2 BS 6.22
3 BT 5.87
4 FM 5.87
5 SI_EK 5.49
6 MU 5.16
7 FB 5.04
8 HO 3.99
9 JM 3.37
10 RA 3.06
write_csv(order, paste0(home, "/output/ranked_temps.csv"))nareas <-length(unique(order$area)) +2# to skip the brightest colors that are hard to readcolors <-colorRampPalette(brewer.pal(name ="RdYlBu", n =10))(nareas)[-c(6,7)]
Plot VBGE fits
# Sample 30 IDs per area and plot their data and VBGE fitsset.seed(4)ids <- IVBG %>%distinct(ID, .keep_all =TRUE) %>%group_by(area) %>%slice_sample(n =30)
Can we fit a single Sharpe Schoolfield with area-specific parameters with brms?
# Again, here are the data we are fitting:ggplot(dat, aes(temp, rate, color =factor(area, order$area))) +geom_point(data = dat, aes(temp, rate, color =factor(area, order$area)), size =0.6)
hist(dat$rate)
# Here's the equation (from the r package rTPC):# > sharpeschoolhigh_1981# function (temp, r_tref, e, eh, th, tref) # {# tref <- 273.15 + tref# k <- 8.62e-05# boltzmann.term <- r_tref * exp(e/k * (1/tref - 1/(temp + # 273.15)))# inactivation.term <- 1/(1 + exp(eh/k * (1/(th + 273.15) - # 1/(temp + 273.15))))# return(boltzmann.term * inactivation.term)# }# Add in fixed parametersdat$bk <-8.62e-05dat$tref <-8+273.15# (better visualization including bounds further down)n =10000hist(rnorm(0.8, 1, n = n), main ="e", xlim =c(0, 5))
hist(rnorm(2, 1, n = n), main ="eh", xlim =c(0, 7))
hist(rnorm(0.25, 0.5, n = n), main ="r_tref", xlim =c(0, 2.5))
hist(rnorm(12, 5, n = n), main ="th", xlim =c(5, 30))
# These work nicely with the pp_checkprior <-c(prior(normal(0.8, 1), nlpar ="e", lb =0),prior(normal(2, 1), nlpar ="eh", lb =0),prior(normal(0.3, 0.5), nlpar ="rtref", lb =0),prior(normal(10, 5), nlpar ="th", lb =0))# Now to a prior predictive checkfit_prior <-brm(bf(rate ~ (rtref *exp(e/bk * (1/tref -1/(temp +273.15)))) / (1+exp(eh/bk * (1/(th +273.15) -1/(temp +273.15)))), rtref + e + eh + th ~1, nl =TRUE),data = dat,cores =2,chains =2,iter =1500,sample_prior ="only",seed =9,prior = prior)
Compiling Stan program...
Start sampling
# Global prior predictive check in relation to data. Doesn't loo too informative...dat %>%data_grid(temp =seq_range(temp, n =51)) %>%ungroup() %>%mutate(bk =8.62e-05,tref =8+273.15) %>%add_epred_draws(fit_prior) %>%ggplot(aes(temp, y = .epred)) +stat_lineribbon(.width =c(0.75), alpha =0.3, color ="black", fill ="black") +geom_point(data = dat, aes(temp, rate, color =factor(area, order$area)), size =0.6) +scale_color_manual(values = colors, name ="Site") +labs(y ="Expectation of the posterior predictive distribution", x ="Temperature (°C)") +NULL
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ggsave(paste0(home, "/figures/supp/prior_predictive_check.pdf" ), width =15, height =11, units ="cm")# Now make sure it converges with real data but no random effectsfit_global <-brm(bf(rate ~ (rtref *exp(e/bk * (1/tref -1/(temp +273.15)))) / (1+exp(eh/bk * (1/(th +273.15) -1/(temp +273.15)))), rtref + e + eh + th ~1,nl =TRUE),data = dat,cores =4,chains =4,iter =4000, seed =9,sample_prior ="yes",family = student,prior = prior )
Compiling Stan program...
Start sampling
plot(fit_global)
pp_check(fit_global)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Ok, seems to work with priors and everything. They are roughly as broad as can be and still have a fit. Now fit the full model, with random area effects, more iterations and chains!
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
Compiling Stan program...
Start sampling
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
Compiling Stan program...
Start sampling
Warning: There were 2 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
Compiling Stan program...
Start sampling
Warning: There were 16 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
Compiling Stan program...
Start sampling
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
ts <- dat %>%data_grid(temp =seq_range(temp, n =100)) %>%mutate(bk =8.62e-05,tref =8+273.15) %>%add_epred_draws(fit_global, re_formula =NA, ndraws =8000)
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
p1 <-ggplot(ts) +geom_line(aes(temp, y = .epred, group = .draw), alpha = .1) +coord_cartesian(ylim =c(0.1, 0.4))p2 <-ggplot(ts) +geom_line(aes(temp, y = .epred, group = .draw), alpha = .1) +coord_cartesian(ylim =c(0.1, 0.4))p1 + p2
# Compute quantiles across spaghettiest_opt <- ts %>%group_by(.draw) %>%filter(.epred ==max(.epred)) %>%ungroup() %>%filter(!temp ==min(temp)) %>%# remove values where there was no optimumfilter(!temp ==max(temp))
# Make the main plot (conditional effect of temperature, with and without random effects)# Predictions without random effectsglob <- dat %>%data_grid(temp =seq_range(temp, n =100)) %>%mutate(bk =8.62e-05,tref =8+273.15) %>%add_epred_draws(fit_global, re_formula =NA) %>%#ungroup()
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
group_by: one grouping variable (area)
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
ggsave(paste0(home, "/figures/sharpe_school_bayes.pdf"), width =17, height =17, units ="cm")#ggsave(paste0(home, "/figures/for-talks/sharpe_school_bayes.pdf"), width = 14, height = 14, units = "cm")
Area-specific T_opts as a ridgeplot
# Calculate T_opt by areatsa <- dat %>%data_grid(area =unique(dat$area),temp =seq_range(temp, n =100)) %>%mutate(bk =8.62e-05,tref =8+273.15) %>%add_epred_draws(fitbs, re_formula =NULL, ndraws =10000)
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
t_opt_a <- tsa %>%group_by(.draw, area) %>%filter(.epred ==max(.epred)) %>%ungroup() %>%filter(!temp ==min(temp)) %>%# remove values where there was no optimumfilter(!temp ==max(temp))
mutate: new variable 'above' (character) with 2 unique values and 0% NA
group_by: one grouping variable (area)
count: now 16 rows and 3 columns, one group variable remaining (area)
pivot_wider: reorganized (above, n) into (N, Y) [was 16x3, now 10x3]
mutate (grouped): new variable 'prop' (double) with 7 unique values and 40% NA
# A tibble: 10 × 4
# Groups: area [10]
area N Y prop
<chr> <int> <int> <dbl>
1 FM 10 27 0.730
2 BT 7 15 0.682
3 SI_EK 14 10 0.417
4 JM 37 23 0.383
5 FB 44 3 0.0638
6 BS 16 1 0.0588
7 HO 29 NA NA
8 MU 18 NA NA
9 RA 14 NA NA
10 SI_HA NA 38 NA
# Without impacted?left_join(dat, t_opt_a_sum, by ="area") %>%filter(!area %in%c("BT", "SI_EK")) %>%filter(temp < median) %>%nrow()
left_join: added 3 columns (median, lwr, upr)
> rows only in x 0
> rows only in y ( 0)
> matched rows 306
> =====
> rows total 306
filter: removed 46 rows (15%), 260 rows remaining
filter: removed 92 rows (35%), 168 rows remaining
[1] 168
dat %>%group_by(area) %>%summarise(n =n()) %>%arrange(desc(n))
group_by: one grouping variable (area)
summarise: now 10 rows and 2 columns, ungrouped
# A tibble: 10 × 2
area n
<chr> <int>
1 JM 60
2 FB 47
3 SI_HA 38
4 FM 37
5 HO 29
6 SI_EK 24
7 BT 22
8 MU 18
9 BS 17
10 RA 14